Code
library(scran)
library(scater)
library(scuttle)
library(ggplot2)
library(pheatmap)
library(patchwork)# add tissue types
spe$TissueType <- ifelse(
spe$TumorType == "ccRCC",
yes = "kid", no = "lun")
# exclude unassigned &
# spots tagged for exclusion
nan <- is.na(spe$anno1) | is.na(spe$anno2)
rmv <- grepl("EXCL", spe$anno1)
sub <- spe[, !(nan | rmv)]
# subset regions of interest
ids <- c("INFL", "TLS", "LN")
sub <- sub[, sub$anno1 %in% ids]
sub$anno1 <- droplevels(sub$anno1)
# subset samples of interest
ids <- c("B04_17776", "B06_24137", "B06_24784")
sub <- sub[, sub$sample_id %in% ids]
sub$sample_id <- droplevels(sub$sample_id)
# simplify annotations
sub$anno3 <- as.character(sub$anno2)
sub$anno3[grep(".*ETLS", sub$anno3)] <- "E_TLS"
sub$anno3[grep(".*MTLS", sub$anno3)] <- "M_TLS"
sub$anno3 <- factor(sub$anno3, exclude = NULL)
table(sub$sample_id, sub$anno3)
E_TLS INFL LN M_TLS
B04_17776 68 140 0 94
B06_24137 235 62 33 166
B06_24784 89 173 0 53
E_TLS INFL LN M_TLS
INFL 1 375 0 5
TLS 391 0 0 308
LN 0 0 33 0
# keep features detected in at least 20 spots in any sample
# and spots with at least 200 detected features overall
idx <- split(seq(ncol(sub)), sub$sample_id)
gs <- sapply(idx, \(.) {
y <- counts(sub[, .]) > 0
rowSums(y) >= 20
})
fil <- sub[rowAnys(gs), ]
fil <- fil[, colSums(counts(fil) > 0) >= 200]
cbind(spe = dim(spe), sub = dim(sub), fil = dim(fil)) spe sub fil
[1,] 16643 16643 8814
[2,] 27595 1113 1065
# split by sample
idx <- split(seq(ncol(sub)), sub$sample_id)
lys <- lapply(idx, \(.) sub[, .])
# feature selection & PCA
lys <- lapply(lys, \(.) {
tbl <- modelGeneVar(.)
hvg <- rowData(.)$hvg <- tbl$bio > 0
runPCA(., subset_row = hvg)
})
# within TLS only
tls <- lapply(lys, \(.) {
. <- .[, .$anno1 == "TLS"]
.$anno3 <- droplevels(.$anno3)
tbl <- modelGeneVar(.)
hvg <- rowData(.)$hvg <- tbl$bio > 0
runPCA(., subset_row = hvg)
})thm <- list(
coord_equal(),
theme(legend.key.size = unit(0.5, "lines")),
guides(col = guide_legend(override.aes = list(alpha = 1, size = 2))))
pal <- hcl.colors(nlevels(sub$anno3), "Set 2")
lapply(lys, plotPCA, color_by = "anno3") |>
wrap_plots() + plot_layout(guides = "collect") &
thm & scale_color_manual(values = pal, drop = FALSE) [1] "IGHG1" "IGKV4.1" "IGHG2" "IGLC1" "IGLV3.1" "IGHG3" "IGHA1"
[8] "IGHG4" "CLU" "IGHJ6" "XBP1" "H4C4"
$B04_17776
$B06_24137
$B06_24784
# top variables by sample
sel <- apply(var, 2, \(.) {
o <- order(., decreasing = TRUE)
names(.)[head(o, 20)]
}, simplify = FALSE)
lapply(names(tls), \(.) {
x <- tls[[.]][sel[[.]], ]
cd <- data.frame(colData(x))
cd <- cd[c("total", "anno3")]
hm <- pheatmap(
main = .,
logcounts(x),
scale = "none",
annotation_col = cd,
show_colnames = FALSE)
print(hm)
})[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=de_CH.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=de_CH.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=de_CH.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Zurich
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] SpatialExperiment_1.10.0 patchwork_1.1.3
[3] pheatmap_1.0.12 scater_1.28.0
[5] ggplot2_3.4.3 scran_1.28.2
[7] scuttle_1.10.2 SingleCellExperiment_1.22.0
[9] SummarizedExperiment_1.30.2 Biobase_2.60.0
[11] GenomicRanges_1.52.0 GenomeInfoDb_1.36.3
[13] IRanges_2.34.1 S4Vectors_0.38.1
[15] BiocGenerics_0.46.0 MatrixGenerics_1.12.3
[17] matrixStats_1.0.0
loaded via a namespace (and not attached):
[1] bitops_1.0-7 gridExtra_2.3
[3] rlang_1.1.1 magrittr_2.0.3
[5] compiler_4.3.2 DelayedMatrixStats_1.22.6
[7] vctrs_0.6.3 pkgconfig_2.0.3
[9] crayon_1.5.2 fastmap_1.1.1
[11] magick_2.7.5 XVector_0.40.0
[13] labeling_0.4.3 utf8_1.2.3
[15] rmarkdown_2.24 ggbeeswarm_0.7.2
[17] xfun_0.40 bluster_1.10.0
[19] zlibbioc_1.46.0 beachmat_2.16.0
[21] jsonlite_1.8.7 rhdf5filters_1.12.1
[23] DelayedArray_0.26.7 Rhdf5lib_1.22.1
[25] BiocParallel_1.34.2 irlba_2.3.5.1
[27] parallel_4.3.2 cluster_2.1.5
[29] R6_2.5.1 RColorBrewer_1.1-3
[31] limma_3.56.2 Rcpp_1.0.11
[33] knitr_1.44 R.utils_2.12.2
[35] Matrix_1.6-1 igraph_1.5.1
[37] tidyselect_1.2.0 rstudioapi_0.15.0
[39] abind_1.4-5 yaml_2.3.7
[41] viridis_0.6.4 codetools_0.2-19
[43] lattice_0.22-5 tibble_3.2.1
[45] withr_2.5.1 evaluate_0.21
[47] pillar_1.9.0 generics_0.1.3
[49] RCurl_1.98-1.12 sparseMatrixStats_1.12.2
[51] munsell_0.5.0 scales_1.2.1
[53] glue_1.6.2 metapod_1.8.0
[55] tools_4.3.2 BiocNeighbors_1.18.0
[57] ScaledMatrix_1.8.1 locfit_1.5-9.8
[59] cowplot_1.1.1 rhdf5_2.44.0
[61] grid_4.3.2 DropletUtils_1.20.0
[63] edgeR_3.42.4 colorspace_2.1-0
[65] GenomeInfoDbData_1.2.10 beeswarm_0.4.0
[67] BiocSingular_1.16.0 HDF5Array_1.28.1
[69] vipor_0.4.5 cli_3.6.1
[71] rsvd_1.0.5 fansi_1.0.4
[73] S4Arrays_1.0.6 viridisLite_0.4.2
[75] dplyr_1.1.3 gtable_0.3.4
[77] R.methodsS3_1.8.2 digest_0.6.33
[79] ggrepel_0.9.3 dqrng_0.3.1
[81] farver_2.1.1 rjson_0.2.21
[83] htmlwidgets_1.6.2 htmltools_0.5.6
[85] R.oo_1.25.0 lifecycle_1.0.3
[87] statmod_1.5.0